{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pickle\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "import cartopy.crs as ccrs\n",
    "import pandas as pd\n",
    "import cartopy.crs as ccrs\n",
    "from pyproj import Geod\n",
    "from scipy import stats\n",
    "from Load_and_Process_Datasets import *\n",
    "import geopandas as gpd\n",
    "from palettable.colorbrewer.diverging import *\n",
    "from palettable.colorbrewer.sequential import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Lat, Lon, and Elevation Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x2ad990785c50>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR0AAAEeCAYAAAC3/sHVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOydZ1hUV9eGnxmKNAVpKgICKqFosCBiQxGDGAtKVOzB/tpbLEQ0sX6xQ4yJoqLYYn0VBRUBGyKIDVBBQZBeJPQ+7Xw/eCEizDAzTIV9X5c/PGeXNcPMM/usvfZaNIqiQCAQCJKCLm0DCARC24KIDoFAkChEdAgEgkQhokMgECQKER0CgSBRiOgQCASJosjrpomJCZWWliYpWwgEQiuDoija19dovOJ0aDQaReJ4CASCMNBotCZFhzxeEQgEiUJEh0AgSBQiOgQCQaIQ0SEQCBKFiA6BQJAoRHQIBIJEIaJDIBAkChEdAoEgUYjoEAgEiUJEh0AgSBQiOgQCQaIQ0SEQCBKFiA6BQJAoRHQIBIJEIaJDIBAkChEdAoEgUYjoEAgEiUJEh0AgSBQiOgQCQaIQ0SEQCBKFiA6BQJAoRHQIBIJEIaJDIBAkChEdAoEgUYjoEAgEiUJEhyBRqqqq8Pfff4NUjm278KxlTmi9ZBVVYcPVWGQUVcJYWx37p3yLzpqqYp2ToigsXboU586dQ3l5OVxcXPDu3Ts4OzuDTie/f20FUsu8DbL37nv8+TC50fWVI3tgrfM3YpmzsrIS27dvx61bt3DmzBk4OTlBWVkZHA4HFy5cgLOzc6M+TCYTdDodCgoKYrGJIF5ILXMCAOB+Ql6TggMAv9//iNySKpHNRVEUXr9+jW3btsHS0hIpKSkIDg5G//79kZeXh7y8PGzYsAHXr19v0K+oqAg7d+5E165dYWZmhm3btiEnJ0dkdhGkC1nptCGqmWz0+uUuWBzubcw7aeDm8qFQUWrZ6iI2NhYrV65ERkYGJk2aBHd3d9jZ2TVql5SUBAcHB2RlZYFOp6O0tBQ2NjZwcHDApk2bUFNTg71796K4uBh+fn44ePAgysrKQFEUrK2tMXDgQPTv35+shmQQbisdIjptiICYLKy6GNPkPYqiQKPVfj50NZRx8scBsDHSEmj8jIwMBAQEICAgAHFxcdi+fTsWLFjQrCD06tULx48fh729PZYuXQoWi4Xjx4/X33/9+jXGjBkDRUVF/PDDD7CwsACHw0FsbCwiIiJQUVGBhQsXYs2aNVBTUxPIZoL44CY6xJHchkgvqGzyOqvsH2QfXwwFDW2omvYDvluC+f7P8WTjyPoVT3Z2NkJDQ6GjowMrKyuYmpoCqPW7nDlzBr6+vvj48SPGjRuHJUuWwNnZGRoaGnzZNWnSJCxfvhzFxcVQVlZGZGRkg/tGRkZQV1fH4cOH8f333zfq//LlS+zZswcDBw7ElStXYGFhIcjbQpAwRHTaENzWrFUpL6FqZgvNoTPx+ZIX2vcbh39ghLtvc0DPeIWDBw8iNjYWo0aNQllZGSIjIxEdHY2SkhJMnz4dpqam2L59O0aOHAklJSWB7Vq0aBE0NTXh7OyM3r1716+46tDV1UVyctN+KADo378/Ll26BD8/PwwdOhTLli3Dhg0boK6uLrAtBPFDHq/aCNVMNob8FoacTx/AKs4FKEDFrB/oSirIv7kXqiZ9oPGtMwpDjqIm5wOUdE2gU5MDLRU6vLy8MGHCBLRr1w4AsGrVKqSlpeHp06c4duwYJk2aJOVX9y8ZGRlYtmwZtLW1cfr0aWmb06Yhj1dtnOB3ucjNykDe3z+jnVEvsMv+QVXKC+iMWQnQaAC99qOg5TAHVamvwakuxyyXGdixYk6jGJqZM2di1KhRuHnzJkaMGCGFV8MdIyMjrFq1Crt27ZK2KQQuENFpI6QXVKI4/Cza27pCa8h0cGoqkX1qBT5f3QZGfiradeoBAKC3U4P6N0Ogq6GMLUtHNhm0Z2dnh/z8/PqVj6yhr6+P7OxsaZtB4AIRnTaCsY4aatLfoKPDjwBqxaXTtF1g5qWApqgMZYN/gwLrdq94bZvLquAAgKWlJXJychAfHw8lJSXo6+tDU1NTpHNUV1fjzp07+P7772X6vZBFiOi0ERx7aoNTVQqF9jr115S0OkNJqzN01JWwaYwlckuqYayjhtHWnVscpyNNFBUVMXHiRDg5OUFdXR35+flwcHDAjz/+iEmTJrUopofNZsPb2xv79+8HjUbD3LlzyaOcgBBHchshISEBo0a7oMuiE/innFF/XdiYHHmitLQUN2/exJEjR1BYWIh169Zh1qxZAsf0fP78GTNnzgSTycThw4ehp6cHGxsbhIaGonfv3mKyXn4hwYFtnB07diAvLw/7D/kg+F0u0gsqW8WqRhAoisKjR49w8OBBREZGwt7eHp07d0bPnj1hY2ODwsJCPH36FFlZWVBVVYWxsTFGjx4NLS0thIeHY+/evZgzZw62bdsGRcXahwQnJyesXbsWY8eOlfKrkz2I6LRhOBwOrKyscOrUKQwaNEja5sgEnz59wps3b5CTk4P3798jLi4OmpqaGDJkCExMTFBdXY3379/jzp07qKiowNChQzFt2jQ4OTnVj8FisdCxY0ekp6ejY8eOUnw1sgkRnTaMn58fjh07hqioqEaBdwThYTKZ0NPTw4cPH9CpUydpmyNzkFPmbZT8/Hx4enri6NGjRHBEjJKSEpydnXH27FmSlEwAiOi0YmpqauDm5oaFCxeib9++0janVfLTTz/hxIkTGDhwID58+CBtc+QCIjqtFA6HgwULFkBfXx/bt2+XtjmtFjs7O8THx8PExAT379+XtjlyAYnTaaV4enoiOTkZoaGhJBWomKHT6VBRUSFpNfiEiE4rJC8vD8eOHUNKSgrXL0IVg43gd7nIKGzZ1nlhOQMHQz7g0z8VMNPTwDpnc2ipKYt0DnmAw2mYGY3NZuPx48cwMzODsbEx8ad9ARGdVsi9e/fg5OQEbW3tJu/HZhRjvv/zBkGC7VUU4ffjAAwwbbpPU1x9kYH11+JQ50ONSC7Auag0jLbqBKuumvCPTEVBGwlENDY2Rnx8PKqqqpCdnQ0PDw8UFBSgqKgIVVVV9Ucl9PT0MGfOHPz4449tdseLrLtbISkpKbC0tGzyXjWT3UhwAKCsmoWpvpGI/lTQqE8Vg40br7NwOCwJATFZqGaykVNchfVX/xWcOigAd+PzcDAksYHgAMA/5QzM93+Oaia7Ra9PFhk6dCguXrwILS0t2NjYwM3NDW/fvkV2djaSk5MRFxeHuLg4HDt2DO/fv8c333yD0aNH4/Tp00hISEB1dbW0X4LEIHE6rZADBw4gMzMThw4danSPV8pSAFBRoiPK0wntFBUQ/C4Xzz8V4mZcNsqqWfVtlBRoYLKF/1z4TOsD1z5dhe4vy1AUBQaD0ewh0IqKCgQGBuLy5ct48+YN0tLS0LlzZ/To0QNWVlawsbGBu7s72rdvLyHLRQ/Jp9OGsLS0xO7du6GkpISNGzdCR+ffQ57cUpbWUc3kwG5XKGg0Gmq4ZHBvieDwY4M8Q6PR+Dp1rq6uDnd3d7i7uwOojW5OT09HUlIS3r17h7Nnz+Lt27fw9vYWt8kSh6x0WinJyclYsWIFBg4ciF9++aX+enMrHUkw2qoTfKb3bbVOZVGQlZWF3r1747fffoOGhgZUVVWhpKQEbW1tDBo0SC4c0+QYRBvk8OHDSExMxOHDh+uvVTPZGLArtMHjkjTQVFHEmfkDW6VTWVScPn0aT548QWVlJSorK8FkMhEXF4dff/0V8+fPr2/3ZSUPWYI8XrVBSkpKGl1TUVKA348DMNU3spETWJKUVLMwz/85Ir6oOEFoiIeHBzw8PBpce//+PYYNG4bMzEyMGDECz58/x/79+6Gurg57e3uMGjUKU6dOlemk9GT3qpXCYDBw9OjRBr+IdQww1calRfZQUZLun7+gnIHgd7lStUHesLCwwN27d1FRUYH169fj2bNnCA4ORmBgIEaNGoUbN27AyMgIGzZsQHl5ubTNbRKRPV5RFIVLly7h3r17KC4uxvnz56GqqioqOwkCcvPmTezbtw/h4eFc2xRXMjB0z32U10hvC3vdd+ZY4dRTavO3RjIzM7F582Y8fPgQvr6+GD16tFTsEPsp89jYWKxZswb9+/cHg8HA3r17RTU0QQiuX7+OKVOm8GyjpaaM8wvsoauhLCGrGmOsQ44OiBpDQ0P4+/vDz88PHh4eOHHihLRNaoDIVjqHDx9GXFwcjh8/jvT0dFhYWKCiokImHVytnfz8fFhaWuLVq1cwNjZutn01k41bsdnwuvGW6za5ONDRUCY+HTGTmJgIZ2dn/PHHHxg3bpxE5xb7SqesrAwqKioAakPCWSwWmEymqIYnCICXlxdmzZrFl+AAtc7lKbZGuLx4EHTUm6/QqapExww7I4yy1BfaRi1VJfg1U3GC0HLMzc1x6NAh/PrrrzKT80dkK52PHz9iyJAhyMzMhJKSEtq3b4/s7Gy5jqiUR2JiYuDi4oL3799DS0vw7ehqJrs+hzIFwP/pJxRU/Pvj8eX5qWomG0P33G90pIIXPfTUsXh4d4y3MSCCIyE4HA6sra1x7NgxODg4SGxeicTp9OvXD0eOHMGgQYNgYGCAJ0+ewMzMTDiLCQJDURQcHR0xffp0LF68WCRjfilCTZ0Ub+rwKC9a8xEIWWbPnj1ITk6Gr6+vxOaUSJyOlZUVkpKSMGjQIAwbNgxBQUFYsWKFKKcg8ODKlSsoKirCggULRDamipICT5GwMdLCk40j64Wps6YKfruT0GB1VIeuhjJGW3cWmW0E/nF1dYWLi4u0zQAgYtExMjJCamoqNmzYgMePH2PJkiWiHJ7Ag7KyMqxbtw4XL15sUTE5YfhamMw7tW+0+uGnaihBfJSWlkJXV1faZgAQsegkJSWhU6dOCAkJQVxcHPT09EQ5PIEHa9asgYuLC4YMGSJtUxqtflp7Ai95ICIigu+NBXEjUtF5/fo1Bg0ahAkTJhDBkSD+/v64f/8+YmNjpW1KPc09lhEkA5vNhpeXFy5cuIAbN25I2xwAYjh7lZ+fjwEDBoh6WEITUBQFX19fbN++HaGhoWSnkNCIwMBABAUF4cWLFzKzEBCp6HTs2BE6OjpISEgQ5bCEJkhOTsaiRYuQn5+P+/fv45tvvpG2SQQZ5NmzZ3Bzc5MZwQFEfODTyckJycnJePPmjSiHJXxFYWEhRo8ejdGjR+PVq1dEcAhcef78ucw9eYg0Tuf9+/fo27cvFBQUUFJSIvFdlLYAk8mEi4sL+vbti/3790vbHIIMU15eDgMDA6SlpUml1rpEygozmUzo6+vD0NAQkZGRohya8D927twJJSUl7NmzR9qmEGSckJAQDBw4UCqCwwuR+nQuXLiA6dOnw8zMDNu2bcO9e/fIgU8REh8fjz///BMxMTFkFUloFk1NzSYTuUkbka50srOzYWxsjLlz56K4uBiTJ09GYWGhKKdo00RHR8PFxQVdu5KtaELzDB06FElJScjLy5O2KQ0Qqeh8//33CAwMhJKSEsLDw2FsbAwrKyvs2LEDSUlJ+PDhA9LT00U5ZZtCRUWFnNwn8I2ysjIsLCzw9u1baZvSAJGKzpgxYxAVFYWsrCyoqKjg0KFDCA0NRVpaGpycnDB+/Hj07t2bfHGEhLxvBEEoKirCu3fvYG9vL21TGiBS0enQoQNmzZqFgwcP1l/r1asXTpw4gfT0dCQmJsLIyEimImfliZcvX6Jv377SNoMgJ4SEhMDBwUHmkrSLPDP3xo0bcfXqVZw8ebLJ+wYGBsjOzhb1tG0CWYy5IMgukZGRGDp0qLTNaITIRadr164ICQnB1q1bcefOnQb3GAwGoqKiZOJQorzB4XDw5s0b9OnTR9qmEOSE2NhY9O7dW9pmNEIsNUjMzc2xZ88eHDt2rMH1goIC0Gg0so0uBHWZALW1taVtCkFOsLOz41kNRFqIrcJnQUEBTE1N8fnz5/rcyQCwdOlS0Gg0HDlyRKhx2yorV66EhoYGdu/eLW1TCHJCXFwcRo4ciQMHDmD27Nmg0yVb50wqZYVHjRqFuXPnYubMmfXXCgsLYWFhgQcPHsDa2lrosdsS5eXlMDIywps3b2BoaChtcwhyRFRUFFavXo2srCwMHjwYAwcOxMyZM9GpUyexzy2RYxBfs3r1ahw8eLBBFnptbW1s2rQJnp6e4py6VREcHIwBAwYQwSEIjL29PSIjIxEWFoZx48bhzJkzuHr1qlRtEqvofP/996ipqUFAQECD60uXLsX79+9x8eJFcU7faggODsb3338vbTMIcgqNRoO5uTlmz56N7t27S2SVwwuRJ/H6EjqdjsOHD8PDwwOjRo2ChoYGgNrI2r///htjxozBwIEDYWpqynUMNpuNV69e4c2bN+jSpQtMTEzwzTffSPz5VJp06tRJ5kLZCfJJUVGRUKWJRInYv7mOjo5wcHDAjh07Glzv378/tmzZguHDh+Pu3buN+mVkZGD58uXQ1dXF3Llz8eDBA3h7e2PcuHHo3LkzZs+eDT8/PyQmJrb6SN2ZM2fixIkTSE1NbXSPoigcOnQIGzZskLxhBLmjc+fOyMnJkaoNYnUk15GTk4PevXsjIiKiUcKp0NBQLFy4EBwOB8bGxtDQ0EBmZiaysrKwYMECrFy5spEvIzU1FcHBwXj06BEiIiKQlZUFBQUFWFpaYvjw4Zg2bRoGDRrUYrtliT179iAkJAQhISH1IQdsNhvu7u7IyMhAVlYWrly50upeN0G0eHp6Ql1dHV5eXmKfi5sjGRRFcf1Xe1s0/Pnnn5SFhQWVm5vb6B6TyaRSUlKoBw8eULdu3aJiYmKosrIyvsfmcDhUZWUlFRkZSe3evZvq1q0bNXLkSOr+/fsUh8MR2WuQJkwmk7Kzs6P+/PPP+mt79uyhHBwcqJqaGurcuXOUhYUFVVBQIEUrCbKOj48PtWzZMonM9T/9aKwrTV2kxCA6FEVR27dvp6ysrKj4+HixigGDwaBOnTpF9ezZkxo8eDB169Ytqrq6WmzzcaOyhkVdf5VJ/R6aSN14nUlVMVgC3f+a+Ph4Sk9PjwoLC6POnDlD6erqUqmpqfX3161bRw0ZMoQqKSkRy+shyD+XL1+m3NzcJDIXN9GRyOPVl+zbtw8+Pj5QVlaGk5MTBg8ejBEjRvB0JgsLm83GlStX4OPjg3fv3mH48OFYtWoVnJycxB4V3VS5XTVlBfxneHcscjDDh9yy+vsURYFGozWoE86NR48eYfz48TA2Nsa5c+caHIvgcDhYtmwZnj17huDgYJlKxk2QDd68eYMpU6bg/fv3Yp9LKsGB3KAoCm/fvsWjR4/w9OlThIaGYurUqdixY4fYUisWFhYiICAA+/btg6qqKtzc3DBu3DjY2NiIfK5qJhtD99znWt9bW10JVYV5SLv9Fxg5SeAwqqDazQYqpn1hPMAZkVvH8ixM9/HjRxgaGjaI9K6D+l898wULFmDWrFkie02E1gGTyYSmpiZycnKgqakp1rlkSnS+prCwEF5eXggMDERgYCC+/fZbsc3F4XAQFhaGoKAgXLt2DSNGjICPj49IzzQFxGRh1cUYnm1KX95CdVosOjrOB02pHarTYlH5IQKswmx4nzyHpa7DhJqbzWajU6dOiImJIcGEhCZxc3ODi4sLFi1aJNZ5pBKRzC/a2tr4888/sW/fPowaNQr3798X21x0Oh3fffcdvL298f79e2hra6NXr1548OCByOZIL6hstg0zPw0qhtZQ6tgFihra0LB2hL6bF9rbTsBGj4kIDAwUau4jR47g22+/JYJD4Iq9vT0uXLggtfllQnTqcHd3x5UrV+Du7o5Hjx6JfT51dXX4+PjgzJkzcHd3R1xcnEjGNdZRa7ZNZVIUFDt2aXS9fR8X/LTvBKbPmQd7Dy9sufEWxZVNP6Z9TVpaGrZv346//vpLYJsJrR8mk4lNmzbh8OHD2Llzp9TskCnRAYDhw4fj0qVLmDJlCvbt24eKigqxzzlq1Cj8/vvvGD9+PPLz81s0VhWDjRomB2rKvKs1aA2bhdKoq6AoToPrGu0UceqjMjpM3okX//XFn3/9hb47QnD1RQbP8TgcDubNm4effvqJFN8jNCIzMxOOjo6IiYnBq1evpJrcS+ZEBwBGjhyJBw8e4Pnz5zAzM8OKFSvw6NEjlJWViW3OadOmYdq0afDw8ACHw2m+QxPEZhRj2N772HAtDpUMNs+2GjbOoCgKFW9C66/pqCujvIYFAFDS7opOM/4PJVFXUBoTjPXX4niueK5du4aysjL89NNPQtlOkD9KSkqwefNm9OjRA1OnTsWRI0dQUFDQoA1FUQgMDIStrS3Gjh2L27dvS31XUyYcyQCQn5+Pq1evIi0tDYWFhejWrRtMTExQXl6O6OhoREVF4dOnT2jfvj1MTU1hYmICQ0NDdO3aFd26dUPPnj3Rs2fPJnd0+IXJZKJHjx4ICAgQOENfcztWTVGT+xH5V3/FllN30KenEaJSCvB3dMMVDbMwC3l//wxtl+VYPHMydkzs1eRYbm5umDBhAjw8PASymyCfUBQFQ0NDjBgxAmvWrEF8fDzu3buH0NBQ7Nq1C5qamggJCcGdO3egpqYGX19fODg4SNRGbo5ksR74FIRbt27B29sbc+bMQbdu3ZCWloZbt26hrKwMFRUVoNPpUFJSQnl5OUpLS5GbmwsWi4WSkhI8fPgQSUlJ9UtId3d3uLu7C1yQTklJCSwWCzo6OgLbH/wuVyDBAYB2nXtAzXI4DmxZi8f3gnD5eeNHKCXtrmjffxxqMt4iJd+lyXEqKioQFhbGNS81ofVBo9Fga2uLIUOGwNbWFra2tpgzZw5evXqFn376CSoqKvjuu++watUqWFpaylS2TpkRHQMDAxgbG2Pz5s0829WV1YiLi8PDhw9x5coVfPvtt/D19YW1tTWCgoJw5MgR/P777/D19RV4+53NZiMmJgZGRkYC9eNnx6opOg73QO6FDRj3H0/MXrAUEckFjdrQVdqDVZQDMz2NJse4d+8eBgwYIHPlYwkth6IonDhxAsnJyaiqqgKdToe6ujqKi4sRGxsLGo2GpUuX1rfv16+fWHd/RYHM+HT69u2LhIQE+Pj48GzXsWNHDB06FEuXLsXly5fx+fNnLF68GDNmzMDSpUsxevRohIeHY/78+XBycsLz588FsuPGjRuYP38+QkJCBOrHz45VU9AUlaA7YSNyHv6NqpyPaOoHiaaoDIrNwDpn80b3kpKSsG7dOsyfP1+o+QmyTUBAAPbu3YsOHTrAxMQERkZGaNeuHUxNTXHt2jVcu3ZN2iYKjMz4dIDa0+NjxoyBk5MT1q5dCzMzM777VlZW4pdffsH58+dx/fp1DBw4ENevX8eqVasQFxcnUA6R8PBwTJ48GREREejRowdffYTx6XxJRfxDlEZexskb97El8D2+fNvL34SgNz0bD25dbtCHzWbDyMgI27Ztw8KFC4WalyC7sNlsmJmZwd/fHyNGjJC2OQIj08GBdZiYmCAiIgIsFgv29vbo27cvduzYgXfv3jXbV01NDfv27cOPP/6IW7duAQAmTZqEvn37IigoSCA7hg0bBi8vLyxevJjvPipKCjj54wDoaigLNFcdapbDQVPtgKAbV/B6y3eYbd8NQ7rrYLZ9N6x3NIZ1N/1Gfeh0OkpKSjBt2jSh5iTINjQaDdnZ2a2uZJNMiQ5QG5189OhR5OTkwMfHBwUFBXBxcYGtrS2OHTuG6upqnv3V1dUb5GR2cXFBaGgojx5Ns3DhQsTExCAzM5PvPjZGWniycST2Tf4W7RQFe2tpNBo62Lnh+sWzSMwrw46JvXB+oT12TOyF8pIi6Os3Fh0ajQYzMzO8fv1aoLkI8gGdToeenp5EDmdKEpkTnToUFBTg4OAAb29vpKamYvfu3bh16xYsLS1x4cIFrrE0urq6yM3Nrf+/kZGRUAF/KioqcHNzw99//y1YPyUFTLE1wuXFg6CjriRQXwXVDqDYLMzxi24Qk5OamsrVsb1u3TqsX79e6Ngigmzzyy+/YOzYsXj58iUqK2s3KyiKAoMh3GO8LCCzovMlCgoKcHZ2RmBgIE6fPg0vLy/8/vvvTbbt06cPYmJ4H7bkl++++w6RkZFC9bUx0kLEJif4TOuDwd353YKnQDGrUcVgYeie+4jNKAZQW7/o60qNiYmJWLp0KY4ePYro6Gg8fPhQKDsJss3ixYvx66+/ws3NDTo6OlBRUYGioiLU1dXRrVs3TJkyBadOnUJOTg6ys7Px8eNH5OXlybQoyYXofIm9vT2qqqrg7Ozc4Lqfnx8uXLiA6urqBr/6NTU1UFYWzs+ioaGBqqoqoW1VUVKAa5+u8PMYwNeqR0nfFHTVDsi/vhvF+TmYcfQxVq5ZBxaLhV69/g0K3Lt3L4YMGQI9PT0cOHAA8fHxculoJPDHvHnzkJaWhsrKShQVFYHBYIDBYCA0NBTjx49HUFAQrK2t0b9/f4wePRq9evWCuro6LCwsZDKhv8zE6fDLq1evoKurCysrqwbXnzx5An9/fxgYGGDt2rX110tLS9GhQweh5urSpQvS09NbZC9QKz5+HnaNknp9DV1JBZ3cd6Aw7Dhyz6xDVlUZUgx64mbAzfpI6+DgYPj4+CA2NhYGBgYtto0gP9BoNKiqqtb/vy4Kf86cOY3aUhSFrVu3YtGiRY1KQEkbmdoy54e6IxLZ2dlo3759/fU1a9ZAQ0MDZWVl2LJlS31U8dGjR/Hq1Sv4+voKPBeLxULHjh2RlpYmknw71Uw2gt/lIr2gEjVMFv54mMKzfd17T6fTcGmRPexMdbBs2TKYmJhg/fr1LbaH0LqpqKiAoaEh4uPj0aVL44wG4kYutsz5QVtbG1OmTMHPP//c4Lqenh6YTCa8vb0bHGPQ1NREaWmpUHMpKirCzs5OaL/O19Q9bn1vQofn+G9RemIu8s6uRVXyiybb02i0//3hUO9cNjMzQ1ZWlkjsIbRu1NXVMWXKFBw5ckTapjRA7kQHAPbv349r167h6tWrqKqqQkVFBTIzM/H58+dGbY2NjSXMLR8AACAASURBVPHmzRsIu2IbPHgwnj592lKTG9C9e3fs27cP4LCwbd0SFN77A4Vhx0GxuNfvqmZyMHTPfVDquiJ55CO0DX755RccO3YM8fHx0jalHrkUHW1tbZw6dQp79uyBrq4uunbtiqysLKxevbpR28GDB4NOpyM4OFioudzc3HD8+HGEh4e31OwGrF69Gtu3b8fRP37HuUv/BcrzkXN2HZgF3OOCymvYOPw0D+9Ts1HN5J06g0AAgK5du2Lnzp2YNWsWampqpG0OADn06XxNVVUVSktLedZnPn36NC5fvozbt2837Muo9bFkFFbCWEcNo607N5kQPSwsDNOmTUNoaKjIE7n/8MMPUFVVxeFjJ2HhthL5D/zRefYBKOuZNNmeVfYPck6tRJ+NF+E3bxDPyhEEAlDrG5w8eTKUlZWxcOFC9O/fX+xJ2YFW5NP5GlVV1WYLwk+dOhVRUVENfCF1CbdWX4rBgZBErLoY0yA25kucnJxgY2ODjAze2fsE5fLly4iOjsbu3buhRDGh//klVM36Q1GrM9c+iu11oajVBZnxLzDr5DO+U5kS2i40Gg0nT56EgYEBtm7dCgMDA3To0AGmpqb4z3/+gzdv3kjUHrkXHX5QU1ODvb09Xr16BaB2F6mp7et/yhmY7/+80aMLi8XC8+fPYW9vLzKbCgoKsGLFCgQEBMDY2BgBAQHooqeN+7cDoKbG+8S6uuUwVMQ/RFk1C/b/F4bnnwpFZhehdaKlpYUDBw7gyZMnKCsrQ0ZGBu7evQsDAwO4uLhg9uzZYDK5+xRFSZsQHQ6HAxaLVR80yCvh1j/lDAS/y21w7ebNm+jduzd0dXVFZtPZs2fx+fNnXLhwAWvXrsXOnTvh7u6OgWa6iPJ0gka72sc8iqLArihGTVYCKj9Gg11VCnWrEahMegZ2dTmqmRxMPRaJJ0kty+1MaDvQ6XRoamrim2++wdatW5GcnIzCwkJMnz5dIsLTqkWHxWJh06ZN6NKlC7KysurTVDSXcOvr+0ePHm2QKEkUzJ07F3///Tc6duwIAwMDbNu2rT7Iq4OKIjy65KH4xk5kHpmN7BP/QWGYL8pe3ETW0fn4fG0H6IrKKAo9BgCgAMw6GY3oT40TgBEIzaGiooL//ve/qK6uxvz584Xe6eUXuXckc4PBYMDV1RUUReGvv/5qULa4uWJ4PtP6wLVP1/r/d+7cGa9evZJIBPDTp0+xZMkSKCsrY/nK1ajs2B2/PytCFbN2lUaxmKjOfIeKt2GoeB+OLnMOQlm/Nu+QihIdUZ5O0FIT7tgHoW1TWVmJkSNH4ttvv8Uvv/yCrl27Nt+JB63WkcyNpKQkJCQkIDAwsFGd9NHWnbnmvdHVUMZo638duRRFoaysrP4oRWlpKXbt2gVLS0toa2uje/fuDU61t4SwsDC4urrCy8sL0dHR+HH2TCwZZ4+LiwahvUrtiRWaohJUTfpAd9w6GK26CCW9f19bXSxPU85wAqE51NTUEBgYiHbt2qF37974z3/+I5bHrVYrOurq6gBqo4q/hlvCLV0NZZz8cUCDbfO0tDR06NAB6urqSExMhIWFBRISEnD27FkkJiZi1qxZcHd3B5stXNxMTU0NAgMD8cMPP8Dd3R3Xrl3DlClTGiTStjHSQvgGR6goNfxz0ZVUGiXcLq9hY97p57j8PAOHw5IQEJNFYnoIfKOrq4vDhw/j06dPyMrKwuTJk4X+bHOj1T1ecTgcXLt2DTt37oShoSHPrIFfnoXiFqdz5swZ7Nq1C8uXL4e3tzd+/vnnBvmIExMTYWtri9zc3GZ3nb7k5cuX+O233xASEgJra2t4eHjA3d2d5+HU558KMfVYJAT9i9SJKYnpIQgCk8mEk5MTnJ2d4eXlJXB/bo9XrU50bty4gUmTJmH27Nk4deqUwGVoviY+Ph5nzpxBeXk5bGxsGuUinj59Onr37t3oLBg3ampqsGTJEty9exebN2/GlClTmswKyI0nSfmYdTJaoNcA1ArPk40jmwx+JBC4ERsbi4kTJ+LTp08C95X5uleiwtXVFaGhofD09ISrqysuXLggdGoLALCyssJvv/3G9b6+vj7fOUsoisKSJUtQWFiIDx8+NDglzy9De+rh8mJ7zPGLRjWT/2yBdaEAXzrICYTmKCoqErgcU3O0OtGh0WhwcnJCREQEVq5ciUGDBiE6OrrexyNqNm/eDCsrK7DZbDg4OGDo0KEwMDBAcXExrl69itLSUpiYmCA5ORk3b95EeXk5njx50iJ77Ex1EOXphKF77qO8hv/nbWFrcxHaLkVFRWjXrp1Ix2y1jmQlJSVoaGhAR0dH5G/al+jr6+PJkycwNTXFuXPn0Lt3b5iYmMDExAR3795FWloa/P39kZKSAk9PT0RGRopEALXUlHF+gb1A1SeErc1FaLuMGjUKL168QE5OjsjGbHU+nToCAwOxZs0aPHv2TCQJuPiFw+EgMTERenp6QpUnFpQvneGdNVXw250EFFQ03uYkPh2CsCxcuBAmJibNVt/9mjbjSAZqfSd2dnbw9PSEm5ubtM2RKLEZxY3OlZHdK0JLePnyJX744QckJycLtDHTpkTn7du3GDduHFJSUkCnt9onSK7wEwpAIAhC9+7dcfPmTVhbW/Pdp83sXgHA8+fPMWTIkDYpOMC/aVEJBFHBZDKhoaEhkrFa1beypKQE27Ztw/r16zFhwgRpm0MgtBpKS0tbVI7pS1rVSmf27NkAgKioqPoT5QQCoeXs2LEDzs7OuH79Ovr169fo+I0gtCrRMTU1haGhIREcAkHErFixAlpaWpg4cSIoioKTkxOcnJzg7u4ucEhKq3Ik//LLLygoKMAff/whcN+m8iVTFPjKoUwgtBUoisKzZ88wZ84cZGVl4fHjx+jfv3+TbVv97lVeXh569eqF8PBwWFhYCNS3qW1mLVUlUKBQUsWqv0a2nmvhN6E9oXXx4cMH7N+/H1evXoWrqyt27NjB84hEq929Kisrwx9//IFDhw5h5cqVAgsOt3zJxVWNA+zqcii35SA7UcUBEeGSL8rLy+Hk5IQFCxbg/fv3zRZD4IXcr3TGjh2L27dvw9PTExs3bhS4tEZzWQSbYrZ9N2wea9nmviTVTDaG7rnfZH5pQSKeSQCj/LFlyxakpKTg/PnzfPdptZkDL168iPPnz+Pdu3cwNjaGq6urQOdEhDkEeTYqrVVm6KtisHHjdRbX5F+CJrRvCkErcRBkg7Nnzwp8DIIbci867du3x4wZMxAQEID09HQoKCjgypUrfPcX9hBka/uS8FMHTNCE9k0hCuEiSJa8vDyUlpbC0tJSJOPJveh8iaamJsaPH4/IyEi++/DKl9wcreVLwu/qozmB5kfARSFcBMlSVVUFZWVlkVWJaFWiw2azceHCBVhZWfHdh1u+ZC1VJWiqNu9nD36XKzerHW6PT/yuPgRJaM8NUQgXQbKYmJhAS0sLr1+/Fsl4cr979SUXL15EUVERPD09BepnY6SFJxtHNjokCQC7ghJwNiqNa9/bb3IR/em+TDhBee0I8XLe8rv6qBNobuPw40SuEy5uzmh+hIsgeczNzZGcnMw1JkcQ5H736ktev36NiRMnIiUlpcW5kevgtWPzJdLOV8NLVL7p3J7nrtNGFwusvxrHdeyv64AVVTBwMOQDUvIrYKangXXO5gLV2iK7V/IFRVHo1KkTXrx4AWNjY777tdo4nS/p27cv9PT0EBAQILI8OnW/7vP8n6OAh/BIIgcxt5VMcz6ZjS4WPB+fAPC9+vhaMCKSC3DnbY5AgsFtZSlKwSZxQKIjPT0dSkpKAgkOL1qV6ADAb7/9hoULF2LMmDFQVVUVyZg2RlrY1MxqABCvE5TX6iC1oIKnqDx4/5nn2Lkl1Xw9NjUnboKs9MSZfoOspERLeno6TExMRDZeq3IkA7U5XW1tbXHgwAGRjptbUt1sG3E5QZv7sid/Luc9QDMHgo111OpXHz7T+mDdd+bwmdYHTzaObPAllYftbhIHJHrKyspEmme81YkOAOzatQs+Pj4oLS0V2ZjNCYpGOwWxOUGb+7IXVfIu/er4jT5fu051q48VTj3h2qdro1WLPGx3y4MwyhtDhw7F8+fPUV7ezI8bn8i16DAYDFy6dAlJSUkNrpubm2P48OG4fPmyyObitV1MowF+Hvzt3ghDc19mbXVlnqIy3saA7zLKvJCH7W55EEZ5o6amBgoKClBRURHJeHLp0yktLcVff/2Fw4cP4/Pnzzh58iR69uzZoI2enh4YDN47Ts3xtTPyz5n9sPT8qwa/pCpKdJyYYws7U9FUfmjKAdrcl9lMT71Zn4wonLfysN0tD8Iob5w+fRrDhw+HoqJo5EIuRefs2bM4evQoAgMDsWLFChgYGDS4/+LFC4SHh2Po0KFCz8HNGbnuu2+w+04CyqprU15UMzlYfSmm3knZkl0TbnP+NbN/s192FSWFZkWlpc5bUcTpiBt5EEZ5IjAwEIcOHUJ4eLjIxpTLOJ3i4mL06NEDz549w19//YWOHTvWH0b77bffcPjwYWzevBkLFy6EkpKSwOPzis2h0YCm3hJdDeUmV0L87po0d4K7JWOLGlmvNkF2r0TH/PnzYWFhgfXr1wvct1XF6WhpacHc3BypqamwtbXFtWvXAADh4eHw9vbGq1evGq1+voTNZuPmzZvQ0NCAqakpzMzMGlSO4OWM5KbB/5QzMO/080ZlfvndTm7OAZpTUi322BZ+kfVqE5KIA2orODo64r///a9Ix5RL0amqqkJsbCzs7e0RHh5ev0u1cuVK/PnnnzwF5/Pnz5g5cyaKiorQoUMHPH/+HH/88Qd+/PHH+jbCOhu51RXnJ3CQHweorH/ZuSGNQD15fa9kDXNzc6SlcT8GJAxyKTqfP3+Gjo4O1NXV0blzZ6SkpCAvLw+fPn3iWnqGw+HA398fnp6emDdvHrZv3w5FRUV0794d/fr1a9BWHM7G5kSltTpAyaOOfJOWlibSwEBATrfMKYqqL4FhY2MDBoOBoKAg9OjRg6uH/datW9i1axeCgoKwe/duKCoqIiMjA+Xl5ejVq1eDts1tjzdFexXe+p2QU9pkYix+5pRXBygJ1JN//vnnH3Ts2FGkY8ql6DCZzHpxodFoGDRoECorK5GSkoKsrKwm+xQXF2Pw4MENTsnGx8ejd+/ejWr4cEt3oauhjH0/fNvkdb8m2n/J7be5TSbG4mdOWdkZEhQSqCf/2NjYICZGsHS+zSGXj1dpaWn1h8+qqqoQEhKCvXv3IioqCnfu3MGCBQsatE9PT4evry9GjBjR6Dq3bPa8nJHjbAyavH7yxwGYdzoaBRXcI4R5OZZbmwOUBOrJP926dUNCQkKDp4uWIpeik5qaCmVlZbx//x4LFiyAs7MzjI2NYWhoiPz8/Pp2WVlZ8PX1xV9//YW1a9c22vYzMjLiujICuDsjeTop+fjD8HIsS9oBymazERwcjBcvXmDr1q0iHZsfPxU5DS7bHDhwAB4eHiITHEBORWf8+PHw9/eHpaUlfHx8sHz5cgCAsrIySktLERgYiBMnTuDx48eYMWMGHj161GR+V0tLS7x8+RKZmZkwNDRskU11/gte6S++RNq/8sHBwTh79izu37+PLl264O3btyIXneYC9Qw0VTFs732ROJmJeIkeDoeDv/76C0+fPhXpuHIZHAjUviEURTVI1jV06FA8e/YMtra2WLhwIaZOnQoNDQ2e4+zbtw8nT57EkydPoKurK7Q9gpay+ToxliR5+PAh3N3dsWPHDjg6OsLMzAxKSkqoqKgQWTqQOrjtXjUV7PjlfUHSZJAdMvGxc+dOXLt2DY8ePUKHDh0E6tvqK3wCwPXr19GzZ89Gu1HN8Z///Adqamo4ePCg0HNv/u8bnI9O57v9ypE9sNSxh8R/jbOystC/f3+cO3cOo0aNqr8+bdo0GBgYYNSoUQgICMChQ4egpiaabfqmIpiD3+XyFGl+RVlUtbgITUNRFFauXIlnz54hMDAQ+vr6fPdtE6IjLNnZ2bCyskJaWprAxfqA2g/+gF2h9eex+EVHXQl+HnYtPrMlCNnZ2Rg2bBhmzZqFX3/9tf5Z/Z9//kHv3r2hp6cHHR0d9OvXT+Q5ib7kcFgSDoQkcr2/7jtzrHDqyfV+Hc2tMKW5omwtUBSFrVu34uLFi4iMjOT7iaBVHYMQNYWFhVBVVRX66H7wu1yBBQcACiqYmO//nOu5qj9n9kd2cZVIhcjAwACRkZEYMWIE+vfvXx9Mqauriw8fPkBDQwOFhYWwtrbG4sWLYW5u3qL5uCGqYEiyQyZ+aDQaduzYgYqKCixfvhwXL15s0XhEdAB4e3tj+fLlQmdHa8kHm9eZLXffyAZnvUTlp9DX18eGDRvg6+vbIIK77pldV1cXEyZMwJ07d8QmOqI6Dd5aI7llkV27dqFPnz44ceJEo7AUQZDL4EBRk5qaioEDBwrdv6UfbG5ntr5+shVlJG9dNjhufP/99zh//jw4HE6L52oKUQVDtsZIbllFVVUVt27dwtatWwWqovs1RHRQu3xksQR/PKpjtHXnZo9BiApRRfLeu3cPzs7OXO+7urqCTqfDz8+vyfsvX75ERUVFi2zgJy9zc7TGSG5ZxtzcHHfu3MHy5ctx584docZo847k5ORkDBw4EPHx8QJ55r/m+adCTD0WCUHfLTVlBVQyBFu58Otk5cXq1auhoaGBnTt3cm0TFhYGT09PREdHA6jdajczM4OxsTFoNBqcnZ0RHBzcIjsEhVcZntYSyS0PREZGYsKECbh27RocHByabEN2r7jg7u6O3r17w8vLq8VjPUnKx6yT0QL10VRVgiIdPI9OfI0odmTevXuHkSNH4tOnT1y3xsvLy9GpUycUFRWBTqfD1NQUlpaWKCwsxMuXL2FkZIQ9e/aAxWLh77//xu3bt1tkU3OQeBzZYsaMGVBSUoK/v3+T98nuFReKiopw9OhRAMDmzZtbFO5ta6INTVUllFTxLyAlVUys/c4cZyJTG3yZeGUoFIWfQk1NDSwWCwwGg6voaGhowMTEBPHx8fj48SMqKysREhICNTU1fPPNN7h06RKcnJygoKAABQUFvH79Gn379m2xbU0hyppbhJZRXV0NT09PvHr1CpGRkQL3b/M+nXv37uH7779HYGBgi8cKfpcrkODUQQMa+TYuLbIXq59i06ZNWLNmDbS0eK8Q7OzsMG/ePMybNw8nTpzAwIEDERQUhFWrVsHGxgaXL1/G7du3sXbtWkycOBF+fn5gs0WfsoKcWJc+dRk3BwwYgMzMTERERAiV9qJNr3QoisKBAwdw9+5dPHv2TOhVTp2f4cqLDKH6G+uoNXnQU1wnzkNCQhAVFYVTp04123bp0qWIiorCjBkzoKOjg0mTJgFA/Yn9kSNHAgD69+8POzs7rFixAiwWC4sWLWqxnV9C4nGkR0lJCY4fP44//vgDnTt3xtatWzF58mShvy9tVnQoioKHhwfevXuH8PBwdOnSRahxmvIzCAKvxyVxnDgvLy/HokWLcOzYMb6OOQwYMAADBgzga2wHBwd4eXnh5MmTIhcdEo8jedLS0uDj4wN/f3+4uLjgypUrfH8WeNFmH698fHyQkJCA8PBwdOvWTagxuPkZ+EWjnYLEt3XXr1+P4cOHw8XFRSzju7i44OnTp0hISBDpuCQeR7IcO3YM/fr1A51Ox+vXr3H+/HmRCA7QhkXnzz//RGlpKXbu3IlPnz4JNQYvPwM/zB9qJtFdl4cPH+L27dvw8fER2xzt27fH77//DgcHB5w/f15k43KLx6kr5xz8LpekPxURFEVh//79uH37Nvbv31+fME9UtFnReffuHY4fPw4GgwE7OztcuHBB4DFa6kcw01NvUX9B2bNnD3799VehDrUKgoeHB8LCwuDl5YU//vhDZON+GUw4084Y7VUUUV7Dxvln6TxTwRIE48WLFwBqNxHEQZuP0wGAmJgYDBgwABUVFVBW5p7n+GsEzaHzJZJOu5CamooBAwYgIyNDZDWp+ZnT0dERa9aswcqVK0U2LklnIV42bdoEOp2O3bt3t2gcbnE6bXal8yXdu3eHioqKQIID8PYz8EIaIfoaGhpgMBgNkp6JGxMTEzx48ADe3t7w9vYW2bhk+1y83Llzh2spJ1FARAdATU0NKIpCSEiIQP24+Rm47SSqKStg3+RvBT5fJCyFhYWorq4GUHty3NDQEG/fvhX7vF9iYmKChw8fwsfHBydOnBDJmGT7XLzk5eUJvbnCD212y/xLdHV1cfPmTcyePRt2dnYYOXIk3Nzc0LVr89vVTVVw6KKpIrW64yUlJfD398fly5cRGRmJPXv24KeffgJQWzmjffv2Yp2/KYyNjXHv3j04ODige/fucHR0bNl4ZPtcrJSWlgqcmlQQiE8Htd766upqVFZW4u7du7h37x4iIiLw6tUrod98SR9ALCsrg5eXF86ePYvRo0djzpw5iImJQW5uLnx8fEBRFJSVlVFSUiKyNKSCEhwcjAULFiAmJgY6OjpCj0N8OuKDw+FAUVERbDa7xRUgiE+HC7Nnz4aamho6dOiAAQMGIDo6GkePHsWoUaOwbNkyocetC+xb4dQTrn26ivVL8PbtW/Tp0wcVFRV4+/Yt/v77b4wZMwY2NjY4cuQI1NTUMGDAAKipqeH9+/dis6M5Ro8eDVdX1xZXnSDpLMRHTU0N2rVrJ9KSM42gKIrrv9rbrZuEhARq0KBB1KRJk6i4uDhq8ODB1Llz56iKigrK0NCQevr0qbRNbJZRo0ZR3t7eTd7jcDhUeXk5df78eUpHR4c6deqUZI37in/++YfS1dWlYmNjWzxWFYNF3XidSf0emkjdeJ1JVTFYIrCwbVNWVkapqqpSNTU1LR7rf/rRWFeauki1IdGhKIoqLy+nVFRUKDabTZ07d44aO3YsRVEUdfr0acre3p7icDhStpA7L168oIyMjPj6kJSWllKVlZUSsIo3586do/T19am7d+9K2xRCEwwePJgKDAxs8TjcRKfNP14BgLq6Ojp06IDc3Fyw2WzQ6bVvy+zZs8FisUQaWStqEhMTMXDgQL62+9u3by/yulbCMHPmTFy5cgWTJ09GWVmZtM0hfMWMGTNw9uxZsY1PROd/GBgYICcnBydOnMCMGTMAAHQ6HYcPH8aGDRuQkSHcCXJxY2JigtTUVGmbITAODg4wMjJCWlqatE0hfMXMmTMREhIits8VEZ3/oauri3HjxqG0tBQTJkyoe7yEvb091q5di3HjxqG8vFzKVjbG0tIS79+/F0sOG3FjYmKClJQUvtpWMdi48ToLh8OSEBCTRc5ZiREtLS0sWLAAhw8fFsv4JE7nf/z888+gKAqOjo6g0Wig0WgICAjAhAkTsG7dOkRGRuL48eNYs2aNtE1tgJaWFvT19fH27VvY2NhI2xyBMDc3R1JSUrPtSJpSycNiscQW00VWOv/D0dERI0eOrIstAACcPn0aQG28wdq1a3H06NH6e7LEypUrMXbsWKFSR0oTS0tLxMXF8WzTXJpSsuIRPRRF4caNG5g4caJYxiei0wT5+fkAgAcPHqCysjakfvDgwWjXrh3CwsKkaVqTrFq1CkePHoWrq2t95QZ5wNHREaGhoTyFnJyzkjwMBgMZGRmwtLQUy/hEdJogNTUV/fr1Q48ePRATU3uKnEajYenSpTh06FCLamSJi3HjxuHkyZOYNGkSzwBAcRXPE4aePXtCWVkZ8fHxXNuQc1aSp127dnytQoWFiE4TFBYWQldXF+rq6mAw/v2VnTVrFqqrq9GnTx+Eh4dL0cKmGT9+PP7v//4Pw4cPx+PHjxvcGz16NBQUFKCoqCgzttNoNPTq1QvJyclc25BzVpKFoiicOnUKWVlZYotKJo7kJmAwGHj27BlKSkoa1DfX0NBAaGgorl+/jh9++AGXLl1q8eFFUTNnzhyUlZXB29u7vghaSUkJnj59irKyMuzbtw8BAQEYNmyYlC2tRVNTE8XF3BNviarmOaFpGAwGHj16hMDAQERFRSExMRGmpqZ48OABevfuLZY5ieg0wbhx4xAREYHKyspGeWFpNBrc3NzQsWNHTJkyBR8/fmy2jIs06Nz53y/jkydPYGdnBzU1Nbi4uIg8abqwUBQFGo2GkpISrm3qzllx270i56xaxqRJk5CXlwc3NzccOHAAFhYW0NXVFeucRHSagE6nw9rammcbR0dH2Nra4v79+3Bzc5OQZfzBZDIb+G6KiorqRUhfX5/nl1xSxMTEwMPDAzU1NVi3bh3PtjZGWghZMxwHQz4gJb8CZnoaWOdsDi01wROoEf4lNTUVz549Q1ZWVoMVvbghotMCxowZg3379sHS0lJsnn5hcHJyqk9nQaPRoKGhUR/YKAtb/gUFBZg0aRK2bNkCDw+P+mMnQG0Q4K3YbDz88BkA4GihDxMddSw5/7J+pRORXIA7b3NInI6QsFgshISE4MCBA5gxY4ZEBQcgotMilixZAhaLheHDh2Pt2rXYtGmTtE0CAPTq1QtlZWXIycmBgYEBTE1N8eHDBwCAsrIyampqpGZbSkoK3NzcMHXqVMybN6/BvdiMYvzoF43iL6qk3n6bCxqAr6WSlBMWjhcvXmDu3LlQVVXFnDlz4OHhIXEbSBIvEZCeno4+ffrg06dPYq+0UFVVhbS0NFRXV8PKyorrQc+ePXsiKCgI5ubmYLPZ0NHRwcuXL6GpqQk9PT2wWCyJ5ksGah/zunfvjm3btmH58uUNdkeqmWwM2XMfBQKW9PGZ1kfkBQlbIxwOB1u2bMGJEydw6NAhTJ8+Xbw5c8A9iRdZ6YgAY2NjTJw4EVOnToWrqyusra1haGgIQ0NDoZauZWVlSEhIgK2tLSiKQlhYGG7cuIFHjx4hOTkZxsbGUFRURFpaGrZs2dJohUVRFBgMRv2jlIKCAnbs2AEHBwdYWlpiM0JfywAACwxJREFU9erVEhccAPWC16dPn0Yf+OB3uQILDkDidPhly5YtCAsLQ1xcHDp16iRVW4joiIhDhw7h/Pnz9dUQs7OzkZOTg06dOsHa2hqLFy/G+PHjG/gvuLFy5UoEBQVBQUEBdDodBgYGcHd3x/z582FjYwNFxdo/m4+PT5Nnl+qSr5ubm9dfW7FiBYyMjODr69vi0iLCQqfTsWLFCqxYsQJjx45Fr169MH36dADCiweJ02kef39/XLx4EVFRUdDT05O2OUR0RIWmpiaWLl3a4BqLxUJaWhqioqKwY8cO/Prrrzh//jysrKy4jvPo0SPcv38fycnJyMvLA4vFgoWFRZNtS0pKUFxcjJcvX8LQ0BD6+vqg0Wi4cuUKpkyZ0mg1MXHiRLGdp+GXRYsWQUNDA+np6fj555/RoUMHjB07VijxIHE6zcNisbBhwwbcu3dPJgQHIKIjVhQVFdG9e3d0794dM2bMwMmTJzF8+HBYWVmhrKwMGhoasLKygouLS70YeHt7w8vLC+3bt2/2lG/fvn3x+PFjLFiwAOnp6WCz2fD29saVK1fqD6vKGsrKyvXOy4EDB2LZsmUYPXo0Rlt3ho6GMtdHrK+dySROhz/CwsJgYmIiUxkIiCNZwqSmpiI1NRXt27dHSUkJLl++jMTERNy/fx8fP36Era0t0tPThapCcf78efzyyy9gMplITU0Vu6NQFAwbNgyrV6/GDz/80OTuFQBoqiri+Bxb5JRUS6y6Rmth7969yMjIEFtuHF5wcySTHMlSZteuXdTQoUMpf39/ytDQkPL19RV6rMTERAoAtXLlShFaKF4uXrxI2dnZUSxWbVL1KgaLuvw8nVpy7gW15OwL6vLzdJJwvQWEh4dT/fv3l8rc4JIjmax0pEx0dDQuXLiAzMxMTJw4EbNmzWrReM7Ozvjtt9/Qr18/EVkoXjgcDkaMGIFJkybJXIK01gCDwYCuri4+ffrUolpjwsBtpUNEhyB1kpOTYWtri0+fPsnkOTZ5Z/LkyRg7dizmzp0r0XlJsT2CzNK9e3cMGTIEd+7ckbYprRJXV1cEBARI24x6iOgQZIIJEyYgKChI2ma0SlJSUmBoaChtM+ohokOQCRwdHfHo0SOZOJDamkhKSsKJEyekcsaKG0R0CDJBjx49wGKxkJ6eLm1TWg1nz57F4MGDsWnTJtja2krbnHpIcCBBJqDRaLC2tsaHDx/QrVs3aZsj11AUhV27dsHPzw8PHjxAr169pG1SA4joEGSGbt26kYqfIuD06dO4fPkyIiIi0KVLF2mb0wjyeEWQGbS1tXnmSybwR3BwMNatWyeTggMQ0SHIEKqqqvV1xgjCUVVVhcePH2Pw4MHSNoUrRHQIMkNeXp7Yk6C1dtauXQsHBwf06NFD2qZwhYgOQWYIDQ2Fk5OTtM2QS6qrq7Fs2TI8fPgQvr6+Mn3Yl4gOQSbIy8tDYWGhzO20yAMfP37EkCFDkJeXh6ioKKEyFEgSIjoEmeD169fo27evTP9CyyIREREYNGgQ5s2bhytXrsjF4ynZMifIBImJiTJVxkdeOHPmDDw9PbFs2TJpm8I3ZKVDkAlKS0vl4lda1oiIiKgvHy0vENEhSJ3q6mqEhobKbFyJrJKUlIT8/Hz06dNH2qYIBBEdglRhsViYNGkSOnXqhCVLlkjbHLni7NmzmDFjRn11EHlBvqwltDqio6ORnp6O2NhYufvySJPPnz/j+PHjCA4OlrYpAkNWOgSpEhERAUdHRyI4AkBRFDw8PDB37lx8++230jZHYMhfmiAVWCwW9u7dC29vb9y4cUPa5sgVmzZtQklJCbZt2yZtU4SCiA5BKpw/fx5Xr17Fy5cvYWRkJG1z5IZDhw7h1q1bCA8Ph5KSkrTNEQoiOgSpkJCQgB9++IEIjgAkJiZi165deP36tcQrO4gS4tMhSIWUlBSYmppK2wy5wtPTEz/99JPcCzVZ6RCkQmZmJoyNjaVthtxw9+5dvH79GmfPnpW2KS2GrHQIUiEzM1OmKhTIChRFobq6utH1devW4eDBg1BTU5OCVaKFiA5B4jCZTOTm5qJr167SNkVmyMvLw4IFC9CtWzd07NgRkyZNQlBQUH11jH79+iE8PFzKVooGIjoEiZOVlYVOnTrJ7e6LqKmsrMSECROgqqqKkJAQ5OXlYcKECfj5558xaNAgPHr0CAcPHsS5c+fw8uVLaZvbYkhZYYLECQoKwsGDBxEWFiZtU6QOg8HAtGnToKqqinPnzjVI7cHhcHDx4kVs3LgR06dPh42NDX7++WeEh4fLhT+M1DInyAyrVq1C586d4enpKW1TpEplZSUmT54MRUVFXL58GSoqKk22++eff+Dh4YHU1FTk5OSgc+fOePv2rcznHuImOmT3iiBxYmJisHXrVmmbIVWqqqowZswYGBsbw8/Pj+ejpq6uLm7evIkHDx6guLiYqzjJC0R0CBJHXV29yR2atgKHw8GsWbNgaGgIf39/0OnNu1bpdHqryR9NHMkEiWNkZIRr166Bw+FI2xSJU1xcDFdXV5SUlMDPz48vwWlttL1XTJA6+/btQ1JSEqZOnYoXL15AEL8hh8PB8+fPERUVhZSUFLDZbDFaKloSEhJgZ2cHExMT3L59G+3atZO2SVKBOJIJUqGiogIHDhyAv78/1NTUsPr/27t7l0byOI7jnxgDImIhkkgUD0EtDKjBNaAYbbQRbAJikcr/YCtrFSxExUqwTKEWIogIEZRIILEwiiAmIGIwiwgqUVACEjHMFcfJHbt7t7u3+zsf3q8yQ4bvhPAm85CZjx8VDAY/O17x8PCgaDSq4+NjnZ2daX19XSUlJSovL9fFxYXcbrcWFxfV0NDwS+Z8enrS6empjo6OlE6nnyP3+PioXC6nfD4vn8+n/v5+VVZWfnU9W1tbCgaDmpyc1PDw8C+Z9aXh7BVeJMuyFIlENDs7q3g8rsbGRtXV1alQKOj+/l6JREKtra3yer2qqalRb2+vWlpa/vxCa25uTmNjYwoEAmpvb1dPT88PBciyLCWTSWUyGV1eXiqZTGpvb0+Hh4eqqqpSc3Oz6uvr5XA4ZLPZ5HA4VFZWJrvdrlgspkgkou3tbbW1tf1tvblcTqOjo1pYWNDy8vKru5/xf0F08OLd3t7q5OREmUxGxcXFKi0tlc/n+8dfENIfz33a2NjQ/v6+wuGwlpaW1NfX97zcsixls1lVVFTIbrd/9v61tTVNTEzo+vpaHo9HLpdLTU1N+vDhg7xe7zfdMH5qakqpVEqhUEiSVCgUtLKyopGREXV3d2tmZkZOp/P7PpBXjujgXYjH4woEAhoaGtLNzY3Oz8+VSqVUKBRUVFSkzs5OBYNBDQ4OyuFw6O7u7nkXbWBg4ItR+hbZbFYej0e1tbXq6OhQOByWy+XS+Pj4mznr9L2IDt6NRCKhWCwml8ul6upqeTweOZ1OXV1dKRqNan5+Xul0WsFgUG63W6FQ6Kf8vSCfz2t3d1c7Ozvy+/3q6ur6CVvzehEd4C8ODg60urqqzc1N+f1+TU9P/98jvTlEB4BRX4sO1+kAMIroADCK6AAwiugAMIroADCK6AAwiugAMIroADCK6AAwiugAMIroADCK6AAwiugAMIroADCK6AAwiugAMIroADCK6AAwiugAMIroADCK6AAwiugAMIroADCK6AAwiugAMIroADCK6AAwiugAMIroADCK6AAwqvhfln+y2Wy/GZkEwFvz6Usv2izLMj0IgHeM3SsARhEdAEYRHQBGER0ARhEdAEb9Ds9bWa6HQ9FxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Load data\n",
    "lat_lon_elevation_path = \"/pl/active/nasa_smb/Data/Density_Profile_Lat_Lon_Elevation.txt\"\n",
    "lat, lon, elevation = np.loadtxt(lat_lon_elevation_path, unpack=True)\n",
    "\n",
    "# Plot Data\n",
    "plt.figure(1)\n",
    "fig1 = plt.figure(figsize=(5, 5))\n",
    "ax = plt.axes(projection=ccrs.SouthPolarStereo())\n",
    "ax.coastlines(resolution='110m')\n",
    "ax.set_extent([-180, 180, -90, -65], ccrs.PlateCarree())\n",
    "plt.scatter(lon, lat, linewidth=2, marker='o', transform=ccrs.Geodetic())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate 1980-2019 Mean MERRA-2 P-E at all sites in meters water equivalent (m.w.e.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = \"/projects/erke2265/SOTC_AIS_SMB/nc_files/\"\n",
    "        \n",
    "# Precipitation\n",
    "sn = xr.open_mfdataset(data_path + \"PRECSN_monthly_*.nc\")\n",
    "ls = xr.open_mfdataset(data_path + \"PRECLS_monthly_*.nc\")\n",
    "cu = xr.open_mfdataset(data_path + \"PRECCU_monthly_*.nc\")\n",
    "\n",
    "# Evaporation \n",
    "evap = xr.open_mfdataset(data_path + \"EVAP_monthly_*.nc\")\n",
    "\n",
    "# Calculate SMB (P-E)\n",
    "smb = (sn['PRECSN'] + ls['PRECLS'] + cu['PRECCU'] - evap['EVAP']) / 12 # convert to yearly average\n",
    "smb = smb.sum(dim = \"time\") * 3.154e+7 / 1000 / (2019 - 1980 + 1) # convert to units of m.w.e/yr\n",
    "\n",
    "# Now get P-E at each lat/lon pair\n",
    "p_e = smb.sel(lon = lon, lat = lat, method = \"nearest\")\n",
    "p_e = p_e.values.diagonal()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-29.379999999999967"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "243.77 - 273.15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Density Profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open a pickle file\n",
    "f = open(\"density_profile_10mfl_ov_rr.pkl\",'rb')\n",
    "density = pickle.load(f)\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot a single profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/projects/erke2265/miniconda/envs/a3d/lib/python3.7/site-packages/ipykernel_launcher.py:8: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  \n",
      "/projects/erke2265/miniconda/envs/a3d/lib/python3.7/site-packages/ipykernel_launcher.py:26: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n"
     ]
    }
   ],
   "source": [
    "for j in range(0, len(lat)):\n",
    "# for j in range(0, 1):\n",
    "    print(j)\n",
    "    # Define index\n",
    "    ind = str(j + 1)\n",
    "\n",
    "    # Plot profile\n",
    "    fig = plt.figure(figsize=(6,10))\n",
    "    plt.plot(density['fdm_density_' + ind], density['fdm_depth_' + ind], 'r', \\\n",
    "             linewidth = 2, label = 'FDM')\n",
    "    plt.plot(density['snowpack_density_' + ind], density['snowpack_depth_' + ind], 'k', \\\n",
    "             linewidth = 2, label = 'SNOWPACK')\n",
    "    plt.plot(density['sumup_density_' + ind], density['sumup_depth_' + ind], 'b', \\\n",
    "             linewidth = 2, label = 'Observations')\n",
    "    plt.grid()\n",
    "    plt.ylim([0, 10])\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.legend()\n",
    "    plt.ylabel(\"Depth [m]\", fontsize = 18)\n",
    "    plt.xlabel(\"Density [kg m$^-$$^3$]\", fontsize = 18)\n",
    "    plt.xlim([200, 650])\n",
    "    plt.savefig('Figures/density_profiles/profile_' + ind + \".png\", format='png', dpi=100)\n",
    "    \n",
    "    # Plot map\n",
    "    plt.figure(1)\n",
    "    fig1 = plt.figure(figsize=(5, 5))\n",
    "    ax = plt.axes(projection=ccrs.SouthPolarStereo())\n",
    "    ax.coastlines(resolution='110m')\n",
    "    ax.set_extent([-180, 180, -90, -65], ccrs.PlateCarree())\n",
    "    plt.scatter(lon[j], lat[j], linewidth=2, marker='o', transform=ccrs.Geodetic())\n",
    "    plt.savefig('Figures/density_profiles/map_' + ind + \".png\", format='png', dpi=100)\n",
    "\n",
    "!zip -r Figures/profiles.zip Figures/density_profiles/\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define Vertical Bins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "depth_interest_low = 0\n",
    "depth_interest_high = 10\n",
    "step = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate depth bin weighted (by layer thickness) density averages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize arrays\n",
    "depth = np.arange(depth_interest_low, depth_interest_high + step, step)\n",
    "SUMup_avg = np.zeros([len(depth) - 1, len(lat)]); SUMup_avg[:] = np.nan\n",
    "SUMup_avg_error = np.zeros([len(depth) - 1, len(lat)]); SUMup_avg[:] = np.nan\n",
    "FDM_avg = np.zeros([len(depth) - 1, len(lat)]); FDM_avg[:] = np.nan\n",
    "SNOWPACK_avg = np.zeros([len(depth) - 1, len(lat)]); SNOWPACK_avg[:] = np.nan\n",
    "\n",
    "# Calculate average in each depth bin\n",
    "\n",
    "# SUMup\n",
    "for k in range(0, len(SUMup_avg)): # Layer\n",
    "    depth_low = depth[k]\n",
    "    depth_high = depth[k+1]\n",
    "    filter_func = np.vectorize(lambda depth: depth >= depth_low and depth <= depth_high)\n",
    "    for j in range(0, len(lat)): # Profile\n",
    "        tmp_density = density[\"sumup_density_\" + str(j+1)]\n",
    "        tmp_depth = density[\"sumup_depth_\" + str(j+1)]\n",
    "        tmp_error = density[\"sumup_error_\" + str(j+1)]\n",
    "        if len(tmp_depth) == 0: # No data at all\n",
    "            SUMup_avg[k,j] = np.nan\n",
    "            SUMup_avg_error[k,j] = np.nan\n",
    "        else: # There is data\n",
    "            tmp_density_filtered = tmp_density[filter_func(tmp_depth)]\n",
    "            tmp_depth_filtered = tmp_depth[filter_func(tmp_depth)]\n",
    "            tmp_error_filtered = tmp_error[filter_func(tmp_depth)]\n",
    "            thickness = np.zeros(len(tmp_density_filtered)); thickness[:] = np.nan\n",
    "            if len(tmp_density_filtered) == 0: # No data in depth bin\n",
    "                SUMup_avg[k,j] = np.nan\n",
    "                SUMup_avg_error[k,j] = np.nan\n",
    "            elif len(tmp_density_filtered) == 1: # Only one data point\n",
    "                SUMup_avg[k,j] = tmp_density_filtered[0]\n",
    "                SUMup_avg_error[k,j] = tmp_error_filtered[0]\n",
    "            else: # There is data in depth bin\n",
    "                for jj in range(0, len(tmp_depth_filtered)):\n",
    "                    if jj == 0:\n",
    "                        thickness[jj] = (tmp_depth_filtered[jj] + tmp_depth_filtered[jj+1]) / 2 - depth_low\n",
    "                    elif jj == len(tmp_depth_filtered) - 1:\n",
    "                        thickness[jj] = depth_high - (tmp_depth_filtered[jj] + tmp_depth_filtered[jj-1]) / 2\n",
    "                    else: \n",
    "                        thickness[jj] = ((tmp_depth_filtered[jj+1] + tmp_depth_filtered[jj]) / 2) - ((tmp_depth_filtered[jj] + tmp_depth_filtered[jj-1]) / 2)\n",
    "                SUMup_avg[k,j] = np.nansum(tmp_density_filtered * thickness / np.nansum(thickness))\n",
    "\n",
    "# SUMup site 96 depth data is not sorted in order, so I throw it out.\n",
    "SUMup_avg[:,95] = np.nan\n",
    "\n",
    "# FDM\n",
    "for k in range(0, len(FDM_avg)): # Layer\n",
    "    depth_low = depth[k]\n",
    "    depth_high = depth[k+1]\n",
    "    filter_func = np.vectorize(lambda depth: depth >= depth_low and depth <= depth_high)\n",
    "    for j in range(0, len(lat)): # Profile\n",
    "        tmp_density = density[\"fdm_density_\" + str(j+1)]\n",
    "        tmp_depth = density[\"fdm_depth_\" + str(j+1)]\n",
    "        if len(tmp_depth) == 0: # No data at all\n",
    "            FDM_avg[k,j] = np.nan\n",
    "        else: # There is data\n",
    "            tmp_density_filtered = tmp_density[filter_func(tmp_depth)]\n",
    "            tmp_depth_filtered = tmp_depth[filter_func(tmp_depth)]\n",
    "            thickness = np.zeros(len(tmp_density_filtered)); thickness[:] = np.nan\n",
    "            if len(tmp_density_filtered) == 0: # No data in depth bin\n",
    "                FDM_avg[k,j] = np.nan\n",
    "            elif len(tmp_density_filtered) == 1: # Only one data point\n",
    "                FDM_avg[k,j] = tmp_density_filtered[0]\n",
    "            else: # There is data in depth bin\n",
    "                for jj in range(0, len(tmp_depth_filtered)):\n",
    "                    if jj == 0:\n",
    "                        thickness[jj] = (tmp_depth_filtered[jj] + tmp_depth_filtered[jj+1]) / 2 - depth_low\n",
    "                    elif jj == len(tmp_depth_filtered) - 1:\n",
    "                        thickness[jj] = depth_high - (tmp_depth_filtered[jj] + tmp_depth_filtered[jj-1]) / 2\n",
    "                    else: \n",
    "                        thickness[jj] = ((tmp_depth_filtered[jj+1] + tmp_depth_filtered[jj]) / 2) - ((tmp_depth_filtered[jj] + tmp_depth_filtered[jj-1]) / 2)\n",
    "                FDM_avg[k,j] = np.nansum(tmp_density_filtered * thickness / np.nansum(thickness))     \n",
    "                \n",
    "# SNOWPACK\n",
    "for k in range(0, len(SNOWPACK_avg)): # Layer\n",
    "    depth_low = depth[k]\n",
    "    depth_high = depth[k+1]\n",
    "    filter_func = np.vectorize(lambda depth: depth >= depth_low and depth <= depth_high)\n",
    "    for j in range(0, len(lat)): # Profile\n",
    "        tmp_density = density[\"snowpack_density_\" + str(j+1)]\n",
    "        tmp_depth = density[\"snowpack_depth_\" + str(j+1)]\n",
    "        if len(tmp_depth) == 0: # No data at all\n",
    "            SNOWPACK_avg[k,j] = np.nan\n",
    "        else: # There is data\n",
    "            tmp_density_filtered = tmp_density[filter_func(tmp_depth)]\n",
    "            tmp_depth_filtered = tmp_depth[filter_func(tmp_depth)]\n",
    "            thickness = np.zeros(len(tmp_density_filtered)); thickness[:] = np.nan\n",
    "            if len(tmp_density_filtered) == 0: # No data in depth bin\n",
    "                SNOWPACK_avg[k,j] = np.nan\n",
    "            elif len(tmp_density_filtered) == 1: # Only one data point\n",
    "                SNOWPACK_avg[k,j] = tmp_density_filtered[0]\n",
    "            else: # There is data in depth bin\n",
    "                for jj in range(0, len(tmp_depth_filtered)):\n",
    "                    if jj == 0:\n",
    "                        thickness[jj] = (tmp_depth_filtered[jj] + tmp_depth_filtered[jj+1]) / 2 - depth_low\n",
    "                    elif jj == len(tmp_depth_filtered) - 1:\n",
    "                        thickness[jj] = depth_high - (tmp_depth_filtered[jj] + tmp_depth_filtered[jj-1]) / 2\n",
    "                    else: \n",
    "                        thickness[jj] = ((tmp_depth_filtered[jj+1] + tmp_depth_filtered[jj]) / 2) - ((tmp_depth_filtered[jj] + tmp_depth_filtered[jj-1]) / 2)\n",
    "                SNOWPACK_avg[k,j] = np.nansum(tmp_density_filtered * thickness / np.nansum(thickness))           \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filter layers that do not have info from both models and observations\n",
    "# Also filter based of elevation of smb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ind = np.zeros(SUMup_avg.shape); ind[:] = np.nan\n",
    "depth = np.arange(step / 2, depth_interest_high, step)\n",
    "\n",
    "# Loop through every profile and layer to calculate filter\n",
    "for k in range(0, len(depth)): # Layer\n",
    "    for j in range(0, len(lat)): # Profile\n",
    "        if ~np.isnan(SUMup_avg[k,j]) and ~np.isnan(FDM_avg[k,j]) and ~np.isnan(SNOWPACK_avg[k,j]):\n",
    "            ind[k,j] = 1\n",
    "\n",
    "# ##### Elevation Filter #####\n",
    "# # Elevation threshold in meters\n",
    "# thresh_high = 5000\n",
    "# thresh_low = 0\n",
    "\n",
    "# plt.hist(elevation)\n",
    "# plt.axvline(x=thresh_high, c='r')\n",
    "# plt.axvline(x=thresh_low, c='r')\n",
    "\n",
    "# for j in range(0, len(lat)): # Loop through all sites\n",
    "#     if elevation[j] > thresh_high or elevation[j] < thresh_low:\n",
    "#         ind[:,j] = np.nan\n",
    "        \n",
    "##### Accumulation Filter #####\n",
    "# Accumulation threshold in m.w.e.\n",
    "# Filter out low accumulation\n",
    "thresh_high = 1\n",
    "thresh_low = .2\n",
    "# qualitative_description = \"High Accumulation\"\n",
    "# Filter out high accumulation\n",
    "thresh_high = .2\n",
    "thresh_low = 0\n",
    "# qualitative_description = \"Low Accumulation\"\n",
    "\n",
    "\n",
    "\n",
    "plt.hist(p_e)\n",
    "plt.axvline(x=thresh_high, c='r')\n",
    "plt.axvline(x=thresh_low, c='r')\n",
    "\n",
    "for j in range(0, len(lat)): # Loop through all sites\n",
    "    if p_e[j] > thresh_high or p_e[j] < thresh_low:\n",
    "        ind[:,j] = np.nan\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Determine which sites are used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sites_used = np.zeros(len(lat)); sites_used[:] = np.nan # Nan means not used, 1 means used\n",
    "for j in range(0, len(lat)):\n",
    "    # If there is non nan in ind[:,j]\n",
    "    if np.nansum(ind[:,j]) > 0:\n",
    "        sites_used[j] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Apply filter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply filter\n",
    "filter_SUMup_avg = SUMup_avg * ind\n",
    "filter_FDM_avg = FDM_avg * ind\n",
    "filter_SNOWPACK_avg = SNOWPACK_avg * ind\n",
    "\n",
    "# Calculate depth averages pth\n",
    "SUMup_mean = np.nanmean(filter_SUMup_avg, axis=1)\n",
    "FDM_mean = np.nanmean(filter_FDM_avg, axis=1)\n",
    "SNOWPACK_mean = np.nanmean(filter_SNOWPACK_avg, axis=1)\n",
    "\n",
    "# Calculate depth standard deviations \n",
    "SUMup_std = np.nanstd(filter_SUMup_avg, axis=1)\n",
    "FDM_std = np.nanstd(filter_FDM_avg, axis=1)\n",
    "SNOWPACK_std = np.nanstd(filter_SNOWPACK_avg, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coastline data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = gpd.read_file(\"/pl/active/nasa_smb/Data/ADD_Coastline_low_res_polygon.shp\")\n",
    "crs_epsg = ccrs.SouthPolarStereo()\n",
    "df_epsg = df.to_crs(epsg='3031')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "# Plot Means\n",
    "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,10))\n",
    "# fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18,10))\n",
    "ax1.plot(SUMup_mean, depth, 'bo--', linewidth = 2, label = 'Observations')\n",
    "ax1.plot(FDM_mean, depth, 'ro--', linewidth = 2, label = 'IMAU-FDM')\n",
    "ax1.plot(SNOWPACK_mean, depth, 'ko--', linewidth = 2, label = 'SNOWPACK')\n",
    "\n",
    "# Plot standard deviations\n",
    "ax1.fill_betweenx(depth, SUMup_mean - SUMup_std, SUMup_mean + SUMup_std, \\\n",
    "                  facecolor = \"blue\", color = \"blue\", alpha = 0.2, label = \"\\u00B1\" + \" \\u03C3\")\n",
    "ax1.fill_betweenx(depth, FDM_mean - FDM_std, FDM_mean + FDM_std, \\\n",
    "                  facecolor = \"red\", color = \"red\", alpha = 0.2, label = \"\\u00B1\" + \" \\u03C3\")\n",
    "ax1.fill_betweenx(depth, SNOWPACK_mean - SNOWPACK_std, SNOWPACK_mean + SNOWPACK_std, \\\n",
    "                  facecolor = \"black\", color = \"black\", alpha = 0.2, label = \"\\u00B1\" + \" \\u03C3\")\n",
    "\n",
    "# Plot Dashed edges\n",
    "ax1.plot(SUMup_mean - SUMup_std, depth, 'b', linewidth = 2)\n",
    "ax1.plot(SUMup_mean + SUMup_std, depth, 'b', linewidth = 2)\n",
    "ax1.plot(FDM_mean - FDM_std, depth, 'r', linewidth = 2)\n",
    "ax1.plot(FDM_mean + FDM_std, depth, 'r', linewidth = 2)\n",
    "ax1.plot(SNOWPACK_mean - SNOWPACK_std, depth, 'k', linewidth = 2)\n",
    "ax1.plot(SNOWPACK_mean + SNOWPACK_std, depth, 'k', linewidth = 2)\n",
    "\n",
    "# Plotting options\n",
    "ax1.grid()\n",
    "ax1.set_xlim([300, 600])\n",
    "ax1.set_ylim([0, 10])\n",
    "ax1.invert_yaxis()\n",
    "ax1.legend(fontsize = 12)\n",
    "ax1.set_ylabel(\"Depth [m]\", fontsize = 18)\n",
    "ax1.set_xlabel(\"Density [kg m$^-$$^3$]\", fontsize = 18)\n",
    "ax1.tick_params(axis=\"x\", labelsize=16)\n",
    "ax1.tick_params(axis=\"y\", labelsize=16)\n",
    "\n",
    "# Plot bias with depth\n",
    "FDM_bias = FDM_mean - SUMup_mean\n",
    "SNOWPACK_bias = SNOWPACK_mean - SUMup_mean\n",
    "\n",
    "ax2.plot(FDM_bias, depth, 'ro--', linewidth=2, label = \"IMAU-FDM\")\n",
    "ax2.plot(SNOWPACK_bias, depth, 'ko--', linewidth=2, label = \"SNOWPACK\")\n",
    "ax2.grid()\n",
    "ax2.set_xlim([-65, 65])\n",
    "ax2.set_ylim([0, 10])\n",
    "ax2.invert_yaxis()\n",
    "ax2.set_ylabel(\"Depth [m]\", fontsize = 18)\n",
    "ax2.set_xlabel(\"Density Bias [kg m$^-$$^3$]\", fontsize = 18)\n",
    "ax2.tick_params(axis=\"x\", labelsize=16)\n",
    "ax2.tick_params(axis=\"y\", labelsize=16)\n",
    "ax2.legend(fontsize = 12)\n",
    "\n",
    "# Plot the number of profiles with depth\n",
    "# Perform count\n",
    "counts = []\n",
    "for j in range(0, len(depth)):\n",
    "    counts.append(np.count_nonzero(~np.isnan(ind[j,:])))\n",
    "\n",
    "# # Plot\n",
    "# ax3.plot(counts, depth, 'ko--', linewidth=2)\n",
    "# ax3.grid()\n",
    "# ax3.set_ylim([0, 10])\n",
    "# ax3.invert_yaxis()\n",
    "# ax3.set_ylabel(\"Depth [m]\", fontsize = 18)\n",
    "# ax3.set_xlabel(\"Number of Profiles\", fontsize = 18)\n",
    "# ax3.set_title(\"Number of profiles that have coincident Obs, FDM, and SNOWPACK with Depth\")\n",
    "plt.savefig(\"Figures/density_profile.pdf\", format='pdf', dpi=100)\n",
    "\n",
    "# Map\n",
    "# Generate figure \n",
    "fig, axs = plt.subplots(1, 1, subplot_kw={'projection': crs_epsg},\n",
    "                        figsize=(10, 10))\n",
    "\n",
    "# Plot coastlines\n",
    "axs.set_extent((-180, 180, -90, -65), ccrs.PlateCarree())\n",
    "axs.add_geometries(df_epsg['geometry'], crs=crs_epsg,\n",
    "                      facecolor='none', edgecolor='black')\n",
    "\n",
    "# Plot SNOWPACK minus Observations\n",
    "plt.scatter(lon*sites_used, lat*sites_used, linewidth=2, c='k', marker='o', transform=ccrs.Geodetic())\n",
    "plt.savefig(\"Figures/density_profile_map.pdf\", format='pdf', dpi=100)\n",
    "print(np.count_nonzero(~np.isnan(lat*sites_used)))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print((SUMup_mean[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print((SUMup_mean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print((SNOWPACK_mean))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print((SNOWPACK_bias))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SNOWPACK_bias.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(FDM_bias)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Investigate standard deviations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SUMup_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SNOWPACK_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FDM_std"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(SNOWPACK_std - SUMup_std, 'b', label='SNOWPACK')\n",
    "plt.plot(FDM_std - SUMup_std, 'r', label='FDM')\n",
    "plt.ylim([-30, 30])\n",
    "plt.grid()\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scatter plot of SUMup vs SNOWPACK in the top meter\n",
    "#### This analysis is different than what I do in `Plot_Top_Meter_Density.ipynb`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "domain = np.array([0, 700])\n",
    "ind_snowpack = []\n",
    "depth_ind = 0 # Ind must be an integer. If ind = x, then you will querry from x to x + 1 m. \n",
    "\n",
    "# Loop through every site\n",
    "for k in range(0, len(lat)):\n",
    "        if ~np.isnan(SUMup_avg[depth_ind,k]) and ~np.isnan(SNOWPACK_avg[depth_ind,k]) and ~np.isnan(FDM_avg[depth_ind,k]):\n",
    "            ind_snowpack.append(k)\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.scatter(SUMup_avg[depth_ind,ind_snowpack], SNOWPACK_avg[depth_ind,ind_snowpack], c='k')\n",
    "# plt.errorbar(SUMup_avg[depth_ind,ind_snowpack], SNOWPACK_avg[depth_ind,ind_snowpack], \\\n",
    "#              xerr=0.1 * SUMup_avg[depth_ind,ind_snowpack], fmt='ko')\n",
    "plt.plot([0, 700], [0, 700], 'k--')\n",
    "plt.xlim([250, 525])\n",
    "plt.ylim([250, 525])\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.grid()\n",
    "plt.xlabel('Observed Surface Density [kg m$^-$$^3$]', fontsize=24)\n",
    "plt.ylabel('SNOWPACK Surface Density [kg m$^-$$^3$]', fontsize=24)\n",
    "\n",
    "slope, intercept, r_value, p_value, std_err = \\\n",
    "    stats.linregress(SUMup_avg[depth_ind,ind_snowpack], SNOWPACK_avg[depth_ind,ind_snowpack])\n",
    "plt.plot(domain, intercept + slope*domain, 'k', label='fitted line')\n",
    "\n",
    "rmse = utilities.calc_rmse(SUMup_avg[depth_ind,ind_snowpack], SNOWPACK_avg[depth_ind,ind_snowpack])\n",
    "bias = utilities.calc_bias(SUMup_avg[depth_ind,ind_snowpack], SNOWPACK_avg[depth_ind,ind_snowpack])\n",
    "\n",
    "print(\"R-squared: %f\" % r_value**2)\n",
    "print(\"p-value: %f\" % p_value)\n",
    "print(\"RMSE: %f\" % rmse)\n",
    "print(\"Bias: %f\" % bias)\n",
    "print(\"Slope: %f\" % slope)\n",
    "plt.savefig(\"Figures/SNOWPACK_top_meter_density_scatter.pdf\", format='pdf', dpi=100)\n",
    "print(len(SNOWPACK_avg[0,ind_snowpack]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate SNOWPACK surface density bias and rmse when observed density is larger than 400 $kgm^{-3}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmse = utilities.calc_rmse(SUMup_avg[0,ind_snowpack][SUMup_avg[0,ind_snowpack] > 400], SNOWPACK_avg[0,ind_snowpack][SUMup_avg[0,ind_snowpack] > 400])\n",
    "bias = utilities.calc_bias(SUMup_avg[0,ind_snowpack][SUMup_avg[0,ind_snowpack] > 400], SNOWPACK_avg[0,ind_snowpack][SUMup_avg[0,ind_snowpack] > 400])\n",
    "print(bias)\n",
    "print(rmse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scatter plot of SUMup vs FDM in the top meter\n",
    "#### This analysis is different than what I do in `Plot_Top_Meter_Density.ipynb`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "domain = np.array([0, 700])\n",
    "ind_fdm = []\n",
    "\n",
    "# Loop through every site\n",
    "for k in range(0, len(lat)):\n",
    "        if ~np.isnan(SUMup_avg[depth_ind,k]) and ~np.isnan(FDM_avg[depth_ind,k]) and ~np.isnan(SNOWPACK_avg[depth_ind,k]):\n",
    "            ind_fdm.append(k)\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.scatter(SUMup_avg[depth_ind,ind_fdm], FDM_avg[depth_ind,ind_fdm], c='red')\n",
    "# plt.errorbar(SUMup_avg[depth_ind,ind_fdm], FDM_avg[depth_ind,ind_fdm], \\\n",
    "#              xerr=0.1 * SUMup_avg[depth_ind,ind_fdm], fmt='ro')\n",
    "plt.plot([0, 700], [0, 700], 'k--')\n",
    "plt.xlim([250, 525])\n",
    "plt.ylim([250, 525])\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.grid()\n",
    "plt.xlabel('Observed Surface Density [kg m$^-$$^3$]', fontsize=24)\n",
    "plt.ylabel('IMAU-FDM Surface Density [kg m$^-$$^3$]', fontsize=24)\n",
    "\n",
    "slope, intercept, r_value, p_value, std_err = \\\n",
    "    stats.linregress(SUMup_avg[depth_ind,ind_fdm], FDM_avg[depth_ind,ind_fdm])\n",
    "plt.plot(domain, intercept + slope*domain, 'r', label='fitted line')\n",
    "\n",
    "rmse = utilities.calc_rmse(SUMup_avg[depth_ind,ind_fdm], FDM_avg[depth_ind,ind_fdm])\n",
    "bias = utilities.calc_bias(SUMup_avg[depth_ind,ind_fdm], FDM_avg[depth_ind,ind_fdm])\n",
    "\n",
    "print(\"R-squared: %f\" % r_value**2)\n",
    "print(\"p-value: %f\" % p_value)\n",
    "print(\"RMSE: %f\" % rmse)\n",
    "print(\"Bias: %f\" % bias)\n",
    "print(\"Slope: %f\" % slope)\n",
    "plt.savefig(\"Figures/FDM_top_meter_density_scatter.pdf\", format='pdf', dpi=100)\n",
    "print(len(FDM_avg[0,ind_fdm]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate IMAU-FDM surface density bias and rmse when observed density is larger than 400 $kgm^{-3}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmse = utilities.calc_rmse(SUMup_avg[0,ind_fdm][SUMup_avg[0,ind_fdm] > 400], FDM_avg[0,ind_fdm][SUMup_avg[0,ind_fdm] > 400])\n",
    "bias = utilities.calc_bias(SUMup_avg[0,ind_fdm][SUMup_avg[0,ind_fdm] > 400], FDM_avg[0,ind_fdm][SUMup_avg[0,ind_fdm] > 400])\n",
    "print(bias)\n",
    "print(rmse)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate colorbar limit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate model differences\n",
    "SNOWPACK_diff = SNOWPACK_avg[0,ind_snowpack] - SUMup_avg[0,ind_snowpack]\n",
    "FDM_diff = FDM_avg[0,ind_fdm] - SUMup_avg[0,ind_fdm]\n",
    "\n",
    "# Retrieve the maximum magnitude difference\n",
    "SNOWPACK_max = np.nanmax(abs(SNOWPACK_diff))\n",
    "FDM_max = np.nanmax(abs(FDM_diff))\n",
    "colorbar_thresh = np.max(np.array([SNOWPACK_max, FDM_max]))\n",
    "print(colorbar_thresh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot SNOWPACK minus Observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate figure \n",
    "fig, axs = plt.subplots(1, 1, subplot_kw={'projection': crs_epsg},\n",
    "                        figsize=(10, 10))\n",
    "\n",
    "# Plot coastlines\n",
    "axs.set_extent((-180, 180, -90, -65), ccrs.PlateCarree())\n",
    "axs.add_geometries(df_epsg['geometry'], crs=crs_epsg,\n",
    "                      facecolor='none', edgecolor='black')\n",
    "\n",
    "# Plot SNOWPACK minus Observations\n",
    "cbar = plt.scatter(lon[ind_snowpack], lat[ind_snowpack], c=SNOWPACK_diff, linewidth=5, \\\n",
    "    marker='o', transform=ccrs.Geodetic(), cmap =  RdBu_5.mpl_colormap)\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label(label='kg m$^-$$^3$', size = 32)\n",
    "cbar.ax.tick_params(labelsize=24)\n",
    "plt.clim(-colorbar_thresh, colorbar_thresh)\n",
    "plt.title(\"SNOWPACK - Observations\", fontsize = 24)\n",
    "plt.savefig(\"Figures/SNOWPACK_top_meter_bias_map.pdf\", format='pdf', dpi=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope, intercept, r_value, p_value, std_err = \\\n",
    "    stats.linregress(SUMup_avg[0,ind_snowpack], SNOWPACK_avg[0,ind_snowpack] - SUMup_avg[0,ind_snowpack])\n",
    "plt.scatter(SUMup_avg[0,ind_snowpack], SNOWPACK_avg[0,ind_snowpack] - SUMup_avg[0,ind_snowpack])\n",
    "plt.scatter(SUMup_avg[0,ind_snowpack], intercept + (SUMup_avg[0,ind_snowpack]) * slope)\n",
    "print(r_value**2)\n",
    "print(p_value)\n",
    "print(slope)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope, intercept, r_value, p_value, std_err = \\\n",
    "    stats.linregress(SUMup_avg[0,ind_fdm], FDM_avg[0,ind_fdm] - SUMup_avg[0,ind_fdm])\n",
    "plt.scatter(SUMup_avg[0,ind_fdm], FDM_avg[0,ind_fdm] - SUMup_avg[0,ind_fdm])\n",
    "plt.scatter(SUMup_avg[0,ind_fdm], intercept + (SUMup_avg[0,ind_fdm]) * slope)\n",
    "print(r_value**2)\n",
    "print(p_value)\n",
    "print(slope)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot FDM minus Observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate figure \n",
    "fig, axs = plt.subplots(1, 1, subplot_kw={'projection': crs_epsg},\n",
    "                        figsize=(10, 10))\n",
    "\n",
    "# Plot coastlines\n",
    "axs.set_extent((-180, 180, -90, -65), ccrs.PlateCarree())\n",
    "axs.add_geometries(df_epsg['geometry'], crs=crs_epsg,\n",
    "                      facecolor='none', edgecolor='black')\n",
    "\n",
    "# Plot SNOWPACK minus Observations\n",
    "cbar = plt.scatter(lon[ind_fdm], lat[ind_fdm], c=FDM_diff, linewidth=5, \\\n",
    "    marker='o', transform=ccrs.Geodetic(), cmap =  RdBu_5.mpl_colormap)\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label(label='kg m$^-$$^3$', size = 32)\n",
    "cbar.ax.tick_params(labelsize=24)\n",
    "plt.clim(-colorbar_thresh, colorbar_thresh)\n",
    "plt.title(\"IMAU-FDM - Observations\", fontsize = 24)\n",
    "plt.savefig(\"Figures/FDM_top_meter_bias_map.pdf\", format='pdf', dpi=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot  Observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate figure \n",
    "fig, axs = plt.subplots(1, 1, subplot_kw={'projection': crs_epsg},\n",
    "                        figsize=(10, 10))\n",
    "\n",
    "# Plot coastlines\n",
    "axs.set_extent((-180, 180, -90, -65), ccrs.PlateCarree())\n",
    "axs.add_geometries(df_epsg['geometry'], crs=crs_epsg,\n",
    "                      facecolor='none', edgecolor='black')\n",
    "\n",
    "# Plot SNOWPACK minus Observations\n",
    "cbar = plt.scatter(lon, lat, c=SUMup_avg[0,:], linewidth=5, \\\n",
    "    marker='o', transform=ccrs.Geodetic(), cmap =  Blues_9.mpl_colormap)\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label(label='kg m$^-$$^3$', size = 24)\n",
    "# plt.clim(-colorbar_thresh, colorbar_thresh)\n",
    "plt.title(\"Observations\", fontsize = 24)\n",
    "plt.savefig(\"Figures/Obs_top_meter_map.pdf\", format='pdf', dpi=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.nanmax(SUMup_avg[0,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.nanmin(SUMup_avg[0,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.nanmean(SUMup_avg[0,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.nanstd(SUMup_avg[0,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope, intercept, r_value, p_value, std_err = \\\n",
    "    stats.linregress(elevation[ind_snowpack], SUMup_avg[0,ind_snowpack])\n",
    "plt.scatter(elevation, SUMup_avg[0,:])\n",
    "plt.plot(elevation, intercept + elevation*slope)\n",
    "print(r_value**2)\n",
    "print(p_value)\n",
    "print(slope * 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope, intercept, r_value, p_value, std_err = \\\n",
    "    stats.linregress(p_e[ind_snowpack], SUMup_avg[0,ind_snowpack])\n",
    "plt.scatter(p_e, SUMup_avg[0,:])\n",
    "plt.plot(p_e, intercept + p_e*slope)\n",
    "print(r_value**2)\n",
    "print(p_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate figure \n",
    "fig, axs = plt.subplots(1, 1, subplot_kw={'projection': crs_epsg},\n",
    "                        figsize=(10, 10))\n",
    "\n",
    "# Plot coastlines\n",
    "axs.set_extent((-180, 180, -90, -65), ccrs.PlateCarree())\n",
    "axs.add_geometries(df_epsg['geometry'], crs=crs_epsg,\n",
    "                      facecolor='none', edgecolor='black')\n",
    "\n",
    "# Plot SNOWPACK minus Observations\n",
    "cbar = plt.scatter(lon, lat, c=SNOWPACK_avg[0,:], linewidth=5, \\\n",
    "    marker='o', transform=ccrs.Geodetic(), cmap =  Blues_9.mpl_colormap)\n",
    "cbar = plt.colorbar()\n",
    "cbar.set_label(label='kg m$^-$$^3$', size = 24)\n",
    "# plt.clim(-colorbar_thresh, colorbar_thresh)\n",
    "plt.title(\"SNOWPACK\", fontsize = 24)\n",
    "plt.savefig(\"Figures/SNOWPACK_top_meter_map.pdf\", format='pdf', dpi=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Check if SNOWPACK bias in surface density (top meter) is lower than that of IMAU-FDM in a statistically significant way"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats.ttest_ind(SNOWPACK_diff, FDM_diff - 1, equal_var=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "a3d",
   "language": "python",
   "name": "a3d"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
